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Finding efficient descriptions of how an environment affects a collection of discrete quantum 
systems would lead to new insights into many areas of modern physics. Markovian, or time-local, 
methods work well for individual systems, but for groups a question arises: does system-bath or 
inter-system coupling dominate the dissipative dynamics? The answer has profound consequences 
for the long-time quantum correlations within the system. We consider two bosonic modes coupled 
to a bath. By comparing an exact solution against different Markovian master equations, we find 
that a smooth crossover of the equations-of-motion between dominant inter-system and system- 
bath coupling exists - but requires a non-secular master equation. We predict a singular behavior of 
the dynamics, and show that the ultimate failure of non-secular equations of motion is essentially a 
failure of the Markov approximation. Our findings support the use of time-local theories throughout 
the crossover between system-bath dominated and inter-system-coupling dominated dynamics. 


I. INTRODUCTION 

A Markovian system is one in which the future time 
evolution depends only on the current state, and not on 
its history [T]. In the context of open quantum systems, 
Markovianity generally implies that the reduced density 
operator obeys a first-order differential equation. This 
class of theory has been developed for many years, is 
applied to a vast range of systems, and provides our un¬ 
derstanding of quantum damping and decoherence [5]. 
Recent work, however, presents it with challenges. The 
development of solid-state quantum emitters, such as sin¬ 
gle and coupled quantum dots Hi] and superconducting- 
qubit cavities Hi, demands theories capable of treat¬ 
ing driven or coupled systems damped by complex struc¬ 
tured baths [THII]. Such theories reveal, among other ef¬ 
fects, the possibility of engineering the reservoirs to con¬ 
trol quantum coherence [Hiia. They show that under 
appropriate conditions both quantum coherence [14] and 
entanglement [HHII] can survive indefinitely, even for 
high temperature baths Hills]. 

These problems do not necessarily elude treatment by 
a time-local theory, i.e., a (Markovian) quantum mas¬ 
ter equation. Such theories accurately reproduce the 
intensity-dependent damping of quantum dots in a struc¬ 
tured reservoir HUOllS], for example, and provide recent 
predictions of bath-induced coherence HI and entangle¬ 
ment [Min]. However, there are several master equa¬ 
tions consistent with, and derivable from, the assumption 
of weak coupling |22[I23| . Furthermore, master equations 
are often postulated phenomenologically, by choice of the 
jump operators in the Lindblad form. For problems with 
multiple oscillators and structured baths this choice is not 
straightforward, with different choices plausible in differ¬ 
ent limits. Nor is it innocent: different forms of master 
equation lead to different behavior [24] . Thus it is impor¬ 
tant to establish which, if any, of the various time-local 
theories is correct. 

In this paper we address this question by studying an 
exactly-solvable model, and comparing the exact solution 


against various time-local theories. We consider a model 
of two bosonic modes, V’a.b, with frequencies cou¬ 
pled to a thermally-occupied bath with spectral density 
J(^). This model has a non-trivial Hamiltonian, multi¬ 
ple degrees-of-freedom, and frequency-dependent damp¬ 
ing, yet is exactly solvable. We consider the general case 
where the natural frequencies coa.b differ and the bath 
couples to a superposition of modes, Paifa + and 

calculate the evolution of the coherence, {iflift)- We 
find a complex behavior with multiple regimes, visible 
in Fig. H reflecting the competing effects of the system 
Hamiltonian and the coupling to the bath. We will com¬ 
pare the exact solution with time-local theories, and so 
identify those which correctly capture such physics. This 
allows us to establish their validity in a generic prob¬ 
lem, and avoids the difficulty inherent in studying only 
approximate theories. 

A physical issue we will address is the appropriate 
form of dissipator for systems with multiple compo¬ 
nents. Two different forms are expected on physical 
grounds [3S]. In the case of the two-oscillator model it 
is clear that at resonance, uia = uJb, the damping can 
depend only on the pattern of coupling to the baths. 
Thus we expect collective decay, described by a Lind¬ 
blad form^Lc = + T*b'4’b] + + Tbi’l], 

where >C[V’] is the standard dissipator with jump oper¬ 
ator if [2]. Far off-resonance, however, we expect indi¬ 
vidual decay terms, L, = 'Ei=a,b^i^bPi\ + 

The first form predicts a non-zero steady-state coher¬ 
ence, while the second predicts this vanishes. We will 
show that neither of these forms is, in general, correct, 
and both make misleading predictions outside of limiting 
cases. Nonetheless, we will find that the general behavior 
can be accurately treated by a time-local theory, specifi¬ 
cally a Bloch-Redfield equation. While one might naively 
have expected some smooth crossover between the limits 
captured by Lc and Li, the real answer is more subtle: a 
smooth interpolation exists for the equations-of-motion, 
but the steady-state is singular at degeneracy. This al¬ 
lows mutually exclusive behavior in different regimes, and 
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implies that some useful effects - specifically the protec¬ 
tion of coherence against the bath - are critically sensitive 
to microscopic parameters. Moreover, while the crossover 
can be treated by a time-local theory, this theory is not 
a Lindblad form with the required positive rates. The 
use of such forms is the subject of ongoing debate, since 
they are not completely positive maps |26j . This means 
that they can lead to unphysical density operators, with 
negative eigenvalues. 

A methodological issue in this debate is the procedure 
of secularization. This amounts to removing from the 
equations-of-motion those terms which are time depen¬ 
dent in the interaction picture. It was used in some of 
the earliest work on quantum damping by Bloch and 
Wangsness m, but was then argued to be unneces¬ 
sary by Redfield [55] as well as Bloch [55]. That po¬ 
sition was challenged by the subsequent Lindblad theo¬ 
rem [30]: as argued by Diimcke and Spohn [23], secu¬ 
larization is required to reach a description where Lind- 
blad’s theorem ensures positivity of the density operator. 
Indeed, Lindblad’s theorem guarantees that the density 
operator will remain positive even when there is entan¬ 
glement with an auxiliary system, a criterion known as 
complete positivity [5]. Secularization, which leads to a 
completely positive theory, is clearly appropriate when 
the interaction-picture time dependence is fast, since off- 
diagonal terms then rapidly average to zero. In our case, 
this is far off-resonance, and secularization indeed leads 
from the Bloch-Redfield equation to the form Li. How¬ 
ever, for a tunable system it may occur that the the time- 
dependence in the interaction picture becomes slow in 
certain regimes, i.e., approaching resonance, so that sec¬ 
ularization becomes inappropriate. An interesting im¬ 
proved version of the secularization procedure is studied 
in Ref. [55] . 

Recently, the necessity of secularization has been ques¬ 
tioned [3ll[35]: Simulations indicate that for time evo¬ 
lution following an initially prepared separable state, 
secularization (and even a Lindblad form for the equa¬ 
tion of motion) can be unnecessary for positivity [33], 
and even complete positivity |34j . in particular for time- 
convolutionless [5] and Nakajima-Zwanzig [HS] [35] ap¬ 
proaches. Stronger statements to this effect have also 
been made by Hell et al. m, noting that a conserva¬ 
tion law [35] is violated by secularized theories — we 
discuss this sum rule in detail further below. The ques¬ 
tion of how the operator form of time-local and non- 
Markovian approaches are related is reviewed by Kar- 
lewski and Marthaler [35]. Our focus in this paper is, 
however, on cases where a time-local description is suffi¬ 
cient. This will enable us to explore the entire parameter 
space of a model, and identify the regions where Bloch- 
Redfield equations predict physical behavior. We will 
show that, although the damping is not of Lindblad form, 
the anticipated unstable behavior does not occur within 
the domain of applicability of the theory - specifically, 
so long as the bath remains Markovian. 

The exact results we present are restricted to only a 


subset of possible initial density matrices. We take as ini¬ 
tial conditions a thermal state of the bath and the ground 
state of the two oscillators. The reduced density matrix 
is then Gaussian at all times, and so completely charac¬ 
terized by its second moments. Thus we will be able to 
establish whether the Bloch-Redfield equations are accu¬ 
rate and physical from the dynamics of those moments 
alone. This does not, however, rule out inaccurate or 
even unphysical behavior for arbitrary (non-Gaussian) 
initial density matrices. 

Within the scope of coupled open quantum systems, 
a particular motivation for our work comes from the 
timely theory of “weak lasing” gni introduced in the 
context of polariton condensates. The idea presented is 
that for modes which are close to resonance the (dissi¬ 
pative) radiative coupling can select which linear com¬ 
bination of modes lases (condenses) first. These works 
started from a phenomenological description of radiative 
coupling, in which collective dissipation terms are intro¬ 
duced by hand. In the following we will see, however, that 
the effects of collective dissipation terms are strongly de¬ 
pendent on whether the individual modes are degenerate 
or not. Our work does not consider the general prob¬ 
lem with both drive and dissipation, but the results we 
present for coupling to a single bath suggest there may be 
a need to re-examine how weak-lasing evolves where ra¬ 
diative coupling selects superpositions of non-degenerate 
modes. 

The remainder of this paper is structured as follows. In 
Sec. [IT] we describe the model. In Sec. [ml we present the 
exact solution, and discuss its behavior. In Sec. |IV| we 
discuss the comparison with the Bloch-Redfield equation 
and the naive Lindblad forms mentioned above. We also 
identify the parameter regimes where the Bloch-Redfield 
equation gives physical behavior. In Sec.|V]we develop an 
alternative to the Bloch-Redfield equation, and show it 
to be an improvement both numerically and analytically. 
In Sec. jVIj we give the generalization of our work to the 
case of multiple baths. Finally, in Sec. VH we give our 
conclusions. 


II. MODEL 

The two bosons and the common bosonic bath are rep¬ 
resented by the Hamiltonian H = Hs + Hsb + Hb- The 
system Hamiltonian is Hs = uJa'>Pl'^a terms 

of bosonic annihilation operators ^i. The bath Hamilto¬ 
nian is Hb = where Ci annihilates a boson in 

mode i. The system-bath coupling takes the form 

Hsb = i‘P*aipi + ‘plipl)Y^giCi + li.c.. ( 1 ) 

i 

The complex coefficients (ft determine which pattern of 
system operators the bath couples to, and gi captures 
the overall coupling to mode i. We will assume the bath 
has a continuous density of states parameterized by the 
spectral density J{v) = ~ 
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Since this model is a linear system of coupled harmonic 
oscillators it is exactly solvable. The exact solution for a 
single harmonic oscillator coupled to a bath m is well- 
known [32] ■ The extension to the case of two identical 
oscillators coupled symmetrically to a bath can be found 
in Ref. [43|, and has been used to test master equation 
approaches m- In this special case the normal modes 
exactly match the pattern of bath coupling. The anti¬ 
symmetric mode then decouples from the bath, imme¬ 
diately reducing the problem to one damped oscillator 
and one undamped one. The dynamics of entanglement 
in this case was studied by Paz and Roncaglia m, who 
showed that the undamped mode allows entanglement to 
persist indefinitely. We consider a more general problem, 
including detuning A = — uJb)j2 ^ 0, which prevents 

such a decoupling and leads to finite lifetimes. 

The existence of finite lifetimes at non-zero A can 
be understood by observing that the model above is 
equivalent to a system of two coupled oscillators, one 
of which is coupled to a bath, i.e., the Hamiltonian 
Hs = -f ujd^d'^d + + H.c. with Hsb = 

'tpl 9A + H.C.. This mapping follows on transforming 
this latter problem to a basis in which Hs is diagonal. 
These two equivalent problems are illustrated schemati¬ 
cally in Fig. We will use the basis of Eq. 0 in the 
following; the results of the other problem can be simply 
extracted by the appropriate rotations. 

In what follows, we consider the time evolution of the 
observables Fijit) = focusing in particular on 

the coherence Fab(t) which, as mentioned earlier, distin¬ 
guishes collective from individual decay. Furthermore, 
these observables fully characterize the density matrix 
for the initial conditions we consider; it is Gaussian, so 
that higher moments are related to the F^j by Wick’s 
theorem. We first present the exact solution and dis¬ 
cuss its observed properties, before considering the (non- 
secularized) Bloch-Redfield (BR) equation of motion. We 
will show both analytically and numerically that this ap¬ 
proach reproduces the exact solution, while either of the 
naive Lindblad master equations fail to reproduce the 
exact results. 



FIG. 1. (Color online) Cartoon of the system we consider: (i) 
two bosonic modes of frequencies uja,udb couple collectively to 
a single bath. As illustrated, we take a super-Ohmic bath 
with an exponential cutoff when an explicit form is required, 
(ii) the equivalent problem of two coupled modes with a bath 
coupling to only one of the modes. 


where Wg = Wf, and vice versa, nsiv) is the Bose-Einstein 
distribution function, and d{C,) = —Aa — C)Ab ~ C) + 
iiF*(C)[|</?aP(wb —C) + |<i5hp(cLia —C)] is ths denominator of 
the retarded Green’s function. Here we have introduced 
Rr(^), the analytic continuation of the damping rate to 
the lower half plane K{() = i f dxJ{x)/{x — C + *0). For 
real C the real part of K{Q is proportional to the spectral 
density, while the imaginary part follows from a Kramers- 
Kronig relation. In the numerical results which follow we 
use the form of spectral density illustrated in Fig. For 
numerical evaluation, it is computationally more efficient 
to write this as a convolution: 


Aj{t) = [ d,T f daDi(t — t)* Dj{t — a)a[a — t), (4) 
do do 

F>i{t) =‘Pi J 


dC (Wi - C)e 


- np-Ft 


27r d{( -h iO) 

q:(t) = J di'J{v)nB(v)e~^''' 


(5) 

( 6 ) 


One may readily check that this is equivalent to 
Eqs. (^). 


III. EXACT SOLUTION 

The exact time evolution can be readily found by us¬ 
ing a Laplace transform to write the system operators in 
terms of the t = 0 bath operators, and then evaluating 
Fij{t) using thermal correlations for the bath operators 
at t = 0. With the oscillators in the ground state at t = 0 
we find: 


Ao{t) = 

J duJ {v)nB {v) W* [v, t)Wj At), 

(2) 


fdC 

^ J 2tt {v — — i0)d{(^ -f iO) ’ 

(3) 


A. Behavior near degeneracy 

Using the above expressions, one may directly find 
how the coherence evolves with time, as the detuning 
A = {(jja — Wb)/2 changes; this is shown in Fig. As is 
clear in this figure, the behavior at degeneracy ^ = 0) 
and away from degeneracy is different. Near degeneracy 
there is strong, long-lived coherence, corresponding to 
the noise-induced coherence recently analyzed for few- 
level models m [T7j. What is not immediately clear 
however is that the degenerate limit is in fact singu¬ 
lar: the form of Fab(po) is discontinuous as a function 
of frequency. We next turn to discuss how and why this 
singular behavior occurs. 























4 



FIG. 2. (Color online) (a) Evolution of coherence Fab with 
time (vertical) and detuning (horizontal), found from the ex¬ 
act solution. We use a super-Ohmic density of states with an 
exponential cutoff, J{uj) = {uijujo)^ ■ This form is 

written such that its peak value is at tj = cuq, and J(a;o) = Jo- 
We choose z = 3, and we measure all energies and times in 
units such that cda = 1. In these units cuq = 0.9, Jo = 0.001, 
and the thermal occupation is controlled by fcsT = 0.52. The 
horizontal axis is the detuning, found by varying ujb for fixed 
LJa- (b) Vertical slice at uJb = 0.9 corresponding to A = 0.05, 
showing comparison between exact and Bloch-Redfield theo¬ 
ries. (c) Horizontal slice at t = 200. In the secularized theory 
the coherence Fab = 0 for all times. 


B. Late time asymptotes 

The long-time asymptotes of the observables can be ob¬ 
tained from the pole structure of the Laplace-transform 
solution |33]. For late times, simplifies signifi¬ 
cantly because the pole at ^ — *0 lies on the real 

axis and so has a vanishing decay rate, while the poles 
of d(C + *0) are generically off axis and so have decayed 
at late times. Thus oo) = — 

v)/d(y + zO). This gives a simplified expression 

F,,(oo) = J dvJ{v)nB{v) ■ (9) 

As can just be seen in Fig. away from the resonance 
point, the off-diagonal coherence decays at late times to 
a small value, but not strictly to zero. However, if the 
bath density of states and occupation are strictly flat, 
i.e. if Jiv) = Jo,nB{v) = no, then one may show that 
the asymptotic value Fij{oo) vanishes for A 0. In this 
case Eq. simplifies considerably, as K(v) = ttJq for 
a flat bath, so div + zO) becomes a simple polynomial. 
This integral then has only four simple poles, and one 
may readily check that it exactly vanishes - except at 
Wq = LOb, where two of the poles coincide and cancel with 
the zeros of the numerator. The small residual coherence 
that exists away from resonance in Fig. [^is thus due to 
the frequency dependence of nB(v), J{v)- 


The origin of the discontinuity at degeneracy is the 
emergence of a slow mode, whose lifetime diverges as 
A —7^ 0. The decay rates of oscillations can be extracted 
from the poles C° of the the Keldysh Green’s function, 
i.e. solutions of d(C°) = 0. Writing Uafi = A leads to 
an expression: 


0 = A2 + A(|:^,P_|^^|2)^^*(^0) 

- (H - C°) - C° - + IPbl^]) . (7) 


At A = 0, the first line vanishes so it is clear that there 
is a pole = H, which is real and so entirely undamped. 
This has a simple physical interpretation: at A = 0, there 
is no coupling between the combination ^^od the 

orthogonal combination of helds. As such, the orthog¬ 
onal combination is entirely undamped, and maintains 
its original state. With non-zero detuning, beating be¬ 
tween the modes ijji=a,b means the orthogonal combina¬ 
tion evolves into Pitpi with time, and is thus damped. 
Considering small A perturbatively gives: 


4A2|v.an<P.P K'{n) 


i\Pa\^ + \Pb\^r\Km^ 


+ OiA^). (8) 


This explains the Gaussian form of the singular response 
visible in Fig. we expect Fab{A,t oo) Fab{A) -I- 
C exp(—aA^t) at large t, where C, a are constant factors 
and Fab{A) is a smooth function. 


IV. BLOCH-REDFIELD APPROACH 

So far, we have seen that the exact solution of the 
bosonic problem does show a crossover between strong 
coherence at degeneracy and weak coherence, due to a 
frequency-dependent spectral density, away from degen¬ 
eracy. However, this crossover occurs as a function of 
time, with coherence surviving over a range at ~ A~^. 
A similar quadratically diverging timescale is found in the 
V-type system M- We now turn to consider whether the 
behavior of the coherence, and other observables, can be 
reproduced by a time-local master equation. 

We may first note that neither of the naive forms (in¬ 
dividual or collective damping) discussed in the intro¬ 
duction reproduce the correct behavior. Separate decay 
predicts that the coherence vanishes, for all times and 
detunings. The collective decay model does predict a 
strong long-lived coherence close to resonance, and in¬ 
deed a quadratically-diverging lifetime. After this time, 
however, the coherence decays to zero, rather than the 
non-zero value predicted by the exact solution. More 
significantly, however, the collective form fails to repro¬ 
duce the behavior once the detuning becomes significant. 
This may be seen from the steady-state populations: 
for large detuning or weak-coupling these correspond 
to equilibrium with the bath, so Faa = nB{oJa),Fbb = 
nB(pJb), whereas the Lindblad form gives equal pop¬ 
ulations. As pointed out by Cresser for the Jaynes- 
Gummings model |45j . such master equations do not 
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reach canonical equilibrium. More generally, since 
is parameterized by one pair of forward/backward rates, 
it cannot account for the presence of two frequencies in 
the dynamics at which the bath should be sampled. As 
can be seen from Fig. this occurs above a critical value 
of the detuning. Thus this model cannot possibly be ac¬ 
curate in this regime, unless the bath and its occupation 
are flat. 

Thus, neither naive form of dissipator can give a full 
account of problems with multiple system frequencies 
and structured baths, particularly if one seeks to ana¬ 
lyze coherence. Unfortunately many interesting prob¬ 
lems in solid-state quantum optics fall in this class, as 
discussed in the introduction. We will now show, how¬ 
ever, that a Bloch-Redfield equation does reproduce the 
correct behavior, as long as one does not secularize the 
final result. Such an approach is frequently stated to be 
invalid, as it leads to negative rates and instabilities. We 
will however show analytically that such instabilities oc¬ 
cur in a much restricted parameter regime, and, in fact, 
only when the Markov approximation breaks down. The 
non-secularized theory is, also, often argued to be in¬ 
valid on the related grounds that it is not a completely 
positive map, and may not even be a positive one. We 
will however show analytically that, although the map is 
not positive, it preserves positivity for almost all Gaus¬ 
sian states. Furthermore, we find numerically that these 
states soon dominate under the time evolution, even if 
dangerous ones are present in the initial conditions. 

Following the standard method [2] one finds the master 
equation has the form: 

dtp = -i[H,p\ + - [p, 

■ij 

+ (2^]M - [p,■ (10) 

b 

Here the Hamiltonian includes Lamb shifts H = Hs — 
J2ij hijPiPj’tpjipj. The matrices h can be written 

in a compact form. 


tions of motion for the quantities derived from these 
master equations. In order to simplify these equations, 
it is convenient to note that the phase of the complex 
coefficients ipi can be eliminated by a phase twist of the 
original operators, and we thus assume pi is real from 
hereon. We may then define the vector of real quantities 
f = and produce an equation of 

motion dti = —Mf -f fo where 


M = 


0 

‘^^a^PbK 

K-2ipa‘PbKa 


0 

2ipb^aKb 


with Fo = (Wf, - ^pIK'^) - {uJa 
PaK'g^ -I- None of these 

bath mode occupations, however 


PaPbK ^aPbK 
V>bPaK -pbPaK 
To —Fq 

Fo Fo 

(13) 

- and Fo = 

rates depend on the 
the constant vector 


fo = plK'b^^, (paPb‘2K[, -(fiaPb^SKi^y in¬ 

volves the excitation rate, so that populations are pro¬ 
portional to the bath occupations as expected. 

The result of time evolving this closed set of equations 
is shown in Fig. ib, c), and clearly compares very well 
to the exact solution. Moreover, we can easily see that 
secularizing this set of equations, as is often claimed to 
be a crucial step [23], could only decrease the agreement: 
secularization can be shown to be equivalent to setting 
all terms involving the product ipaPb to zero, thus remov¬ 
ing the off-diagonal blocks of Eq. (13) and the last two 
elements of the vector fp. This then makes the coherence 
Fab{t) identically zero. This is as expected for a secular 
theory: a non-zero detuning LOa bJb means the master 
equation contains no cross terms between modes a, b and 
thus no coherence arises. Note that the coherence in the 
secular theory is identically zero, whereas that in the ex¬ 
act result decays to a small value after the time 1 /(q;A^), 
see Eq. Q. We can thus identify this timescale as that 
controlling the secular approximation. 


A. Stability of time evolution 




h = 


Ka 

KT^SK 

K±i5K';\ 

) 

(11) 

K 

F" - i6K'\ 

(12) 

K" + iSK' 

j’ 


with the upper (lower) signs in Eq. (11) for (E^). 
Here we have introduced several new pieces of no¬ 
tation. We have used the shorthand Kt = K{uji) 
in terms of the Hilbert transform (analytic continua¬ 
tion) defined previously, and have also defined Hilbert 
transforms of the excitation (absorption) rate = 
i J (i^nB(^) J(^)/(^—Wi-l-jO), and de-excitation (emission) 
rate i J (i^(nB(^) +1) J(^)/(^ —-|-j 0). Note that 

this means Ki = — Primes signify real and imag¬ 

inary parts and X = {Xa + Xb)j2, SX = (Xa — Ff,)/2. 

While Eqs. ([l0}|^ fully describe the equations of mo¬ 
tion, it is more convenient to use the (closed) set of equa- 


The frequently stated reason [23] for secularizing the 
equation of motion is that it is required to ensure the 
equation is of Lindblad form with positive rates, i.e. 
that the master equation take the form p = —i[H,p] + 
Ai(2AjpA|’ - [AjAj,p] + ) with Ai > 0. This is desired 
so that Lindblad’s theorem can guarantee complete pos¬ 
itivity of the density matrix. In addition, negative de¬ 
cay rates may lead to exponentially growing observables. 
Despite its near-perfect match to the exact solution, our 
non-secularized equation clearly fails these requirements. 
Eq. (10) can be put into Li ndb lad form by diagonalizing 
the matrices Eq_ (11), however the eigenvalues 

are X„=K± S, where {K',f + \6Ky^ > (F;)^. 
This means that except when = 0, one rate is al¬ 
ways negative. Despite this, there have been several re¬ 
cent works |32j which suggest it is not established that 
this formal problem leads to any practical difficulties in 
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applying such a theory. 

In our problem, we are able to find precise conditions 
under which the negative rates in the Lindblad form 
cause a practical problem. Operationally, our problem 
is to solve the four linear coupled equation for the com¬ 
ponents Fij. This method will fail if the matrix M has 
negative eigenvalues. For a Gaussian problem such as 
the one we consider here this condition is in fact the only 
practical consideration; all higher moments factorize by 
Wick’s theorem and so positivity of the eigenvalues of M 
ensures the dynamics remains bounded. Remarkably, the 
eigenvalues of M can be found in closed form. They are 
M(/)i = jjLitpi with 


/r* = ^[Ka + iFfe] ± -v/3fJ[Q] ± IQI, 

. . 1 r -|2 

Q = 2KaKh + - ^Ka - Kh+ i{uJa - Wb) 


(14) 


where Ki = ^^Ki. It is clear that when oja = uJb, one 
finds Q = [Ka + which means 3?[Q] -I- |Q| = 

(3fi[Rra -I- Kb])'^. Thus the Bloch-Redfield form recovers 
the fact there is a zero eigenvalue at degeneracy. 

From this closed form we may check that the eigenval¬ 
ues remain positive (stable) as long as 

+ A(< + K'b){Ka^b - K'bK) > 0. (15) 


The first term is always positive, and thus instability re¬ 
quire two conditions: Firstly, it requires that A{K'g^K'^ — 
K^K”) < 0, placing a constraint on the frequency de¬ 
pendence of J{i') — typically an instability is hard to 
achieve if J{v) has only a single peak, but is possible 
for a multi-peaked structure. Secondly, and more im¬ 
portantly, in order for the second term in Eq. (15) to 
dominate, dK{uj)/dio must be large enough — this cor¬ 
responds directly to requiring that the spectral density 
should vary significantly on a scale J{uj), i.e. that the 
memory time of the bath is comparable to the damp¬ 
ing timescale. If such a condition is satisfied, then the 
Markov approximation is a priori invalid. 

To summarize, as long as the Markov approximation is 
valid a priori - i.e. the bath memory time is short com¬ 
pared to damping time - then the eigenvalues of M are 
positive and the solution is stable. This result shows that 
Markovianity is, for this problem, a sufficient condition 
for stability. This is despite the Lindblad matrices 
always having negative eigenvalues, except at resonance. 


B. Comparison to exact solution near degeneracy 


We have already seen the numerical agreement between 
this BR treatment and the exact result in Fig. We 
may note that near resonance one can compare the per¬ 
turbative solution of the exact problem to a perturbative 
expansion of the BR eigenvalues. Starting from Eq. (141, 


and expanding up to quadratic order in A and 5K, one 


finds: 


Mo = 


SifliflAiK'SK" - 5K'K") 


0{A^). (16) 


Recall that 5K depends on the detuning, vanishing at 
least linearly as A —>■ 0, so that the second term is at least 
second-order in A. This eigenvalue can be compared to 
the exact perturbative result by referring back to Eq. 
and noting that = — 23[C°]- The factor of two 

appearing here is because /i corresponds to the eigenvalue 
of the population equation, whereas the pole in Eq. 
gives the decay of fields ipi. 

Comparing Eq. ® to Eq. (16) one sees that the 


leading-order term in K{uj) is correct, but the second 
term in Eq. (16) is not there in the exact solution. The 


second term is however dependent on the derivative of the 
function K. Thus one finds again that the BR theory is 
correct as long as the Markovian approximation holds, 
i.e. as long as the derivative of the density of states is 
sufficiently small. 


C. Positivity of time evolution 


As we have seen, the Bloch-Redfield time evolution 
is stable, and has the correct steady-state, so long as 
the Markovian approximation is justified. This rules out 
the most dramatic pathologies that could arise from the 
negative rates, and suggests that the dynamics will not 
stray far from the correct behavior. This is consistent 
with the essentially perfect agreement seen numerically. 
We now consider a related issue, of the extent to which 
the negative rates lead to unphysical density matrices 
with negative eigenvalues. 

We first summarize some standard definitions [5]. An 
operator is positive if all its expectation values are posi¬ 
tive, and a map is positive if it is between positive opera¬ 
tors. Since density operators are positive the exact time- 
evolution superoperator, which is a map between density 
matrices, is positive. The secularized master equation 
in fact satisfies the stronger criterion of complete posi¬ 
tivity, which corresponds to positivity in the presence of 
arbitrary entanglement with an auxiliary system. 

The map given by Eq. (10) can be shown to be non¬ 


positive specifically because of the negative eigenvalues 
of the Kossakowski matrices To demonstrate this 

we suppose that has a negative eigenvalue, and work 
in its diagonal basis. We denote the field operator corre¬ 
sponding to the unstable (stable) eigenvector by (t/’d), 
so that there will be terms in Eq. ([l0|) of the form 


- [P, 


cj + 


(17) 


with r < 0. In general neither hij nor k will be diagonal 
in this eigenbasis of L^, so that H contains terms 'iplfpj, 
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for a ll pairs i, j S c,d. Similarly, we will have terms in Eq. 
(101 from of the form pipi — [p, +, for all such 

pairs. However, as positivity requires that all positive 
operators are mapped to positive operators, showing it is 
violated only requires us to construct a single counterex¬ 
ample of a positive operator mapped to a non-positive 
operator, and it is possible to do this despite the non¬ 
diagonal nature of these other terms. To construct this 
counterexample we suppose p describes a pure Fock state 
in the diagonal basis of p = |n, TO)(n, to|. This is a 
positive operator, which is mapped by the first term in 
Eq. (171 to 2rn|n — l,m)(n — l,m|. Furthermore, we 


see that no other term in the infinitesimal time-evolution 
superoperator <i)(p) generates this operator. The Hamil¬ 
tonian and anticommutator terms in Eq. (10) conserve 


the total excitation number, while the jump terms from 
increase it. Thus (n — 1, m|$(p)|n — 1, m) = 2rn < 0. 
Since a positive operator X obeys V|'0) : {ilj\X\'tjj) > 0, 
this fact proves that the map has taken a positive op¬ 
erator to a non-positive operator. This proves that the 
map is not positive. It follows immediately that it is not 
completely positive. An analogous argument applies for 
a negative rate in L^. 


While the Bloch-Redfield Eq. (101 is not positive, it 


nonetheless agrees well with the exact solution. This sug¬ 
gests that the operators which are mapped out of the 
physical space, such as the one constructed above, are 
absent from, or at least a negligible contribution to, the 
density matrix. To investigate this, and explore the do¬ 
main of validity of the theory more generally, we consider 
whether the dynamics is positive for Gaussian states. We 
consider specifically the subset of Gaussian states rele¬ 
vant to the dynamics above, where the baths and initial 
conditions are such that Gij = {ipi'i’j) = 0. 


For Gaussian states the density matrix is positive if the 
uncertainty principle is satisfied |46) , which here is equiv¬ 
alent to Fij being positive semi-definite. This follows on 
noting that for two oscillators any normalized linear com¬ 
bination of the operators ipa, 4’b is a lowering operator fj^ 
with corresponding quadratures x = {fj + fi^)/'/2,p = 
—i{fj — fj^)jy/2, and requiring AxAp >1/2 for all such 
quadratures. Positivity of the density matrix can thus 
be checked numerically by calculating the smallest eigen¬ 
value of Fij. In the Bloch-Redheld solution correspond¬ 
ing to Fig. da) we find that there is a brief transient 
period, up to t « 1, where the state violates positivity 
by a tiny amount. Specifically, the smallest eigenvalue of 
Fij reaches —10“^ in this regime, after which it it 

is always positive or zero, with typical values Am "^0.1. 
More generally, the Bloch-Redfield Am agrees with the 
exact result to four decimal places. The error is hardly 
noticeable, except in that it takes the results slightly out¬ 
side the physical regime at early times. 


The behavior discussed above can be understood by 
deriving the condition under which the Bloch-Redfield 
Eq. (101 preserves positivity for Gaussian states. For a 
time increment At the Bloch-Redfield Eq. (101 implies 


a shift in the Fij, Fij —)■ E/® -(- AtRij. Since the time 
evolution of the density matrix is continuous it can only 
become unphysical if F^^'^ has a zero eigenvalue, which 
becomes negative under the perturbation Rij. Such an 
must be of the form 


na ^/UaUbF'*’ 

rib 


(18) 


with nayUb > 0. From the forms of M and fo we cal¬ 
culate the shift matrix elements Rij for the state Fij\ 
We can then calculate the shift in the zero eigenvalue 
perturbatively, and find it to be negative when 


- 2ipa‘Pb^rianb[Kj cos{(j)) - SK^ sin((/)] < 0. (19) 

This condition gives a range of — rib and (p for which 
the minimum-uncertainty Gaussian state, Eq. (18), is 


mapped out of the space of physical states. If the rates 
are not too different, i.e., the Markov approximation is 
well satisfied, then this range is small. 

In summary, the two-mode Bloch-Redheld equation is 
positive for most Gaussian states. The exceptions are 
rare, being the subset of minimum-uncertainty states de¬ 
fined by Eq. (19). Since the dissipation drives the system 


towards safe Gaussian states these dominate the dynam¬ 
ics, even if the others are present in the initial conditions. 
Indeed, a positivity-violating state is present in the ini¬ 
tial condition for Fig. but its effects are transient and 
quantitatively small. 


V. BEYOND THE BLOCH-REDFIELD 
EQUATION 


From the above we may conclude that over a wide 
range of parameters the Bloch-Redfield theory with¬ 
out secularization accurately matches the exact solution, 
while secularization reduces the accuracy. This however 
leaves open an alternate question: does the BR mas¬ 
ter equation, and the corresponding coupled equations of 
motion for Fij{t), represent the best possible time-local 
theory of this problem? In this section we show that a 
better set of time-local equations exists, and involves a 
minor change to the form of the matrix M that appears 
in Eq. (13). 


A. Sum rule violation 


There are two motivations to suggest that an improved 
equation is possible. The first is that, as noted above, 
the BR prediction for the slowest decay rate of coher¬ 
ence, Eq. (16), does not match the rate derived from the 
poles of the exact solution, Eq. (§. The second rea¬ 
son concerns sum rules as discussed in [ 571155 | . These 















state that for operators which commute with the system- 
bath coupling the time evolution of such operators in 
the full dynamics should be equal to that in the ab¬ 
sence of system-bath coupling. For our model, the op¬ 
erators X = — 'fa4’b and obviously commute 

with Eq. Q. As such, their time derivatives should be 
the same as that following from Hs alone. In terms of 
population equations this corresponds to the statement 
that 

/ = (AtA) = l^bl'^Faa + l^al'^Fbb - 


should obey dtl = ^[2i{uja — uJb)(Pa<fbFab]- In the case 
that (fi are real, this means that one should have: 


/ 

V 



T 

M = (Wa 



( 20 ) 


One may however immediately see this does not hold for 
the solution Eq. (13) of our time-local master equation, 
unless Ka = Kb- We next find an alternative time-local 
equation of motion for the observables Fij(t) that both 
satisfies this sum rule, and gives the exact eigenvalues 
near degeneracy. 


B. Schrodinger picture Bloch-Redfleld equation 


The basis of the alternate approach is to consider the 
Born approximation for the equation of motion, before 
making any Markov approximation. We therefore first 
recall the form of the integro-differential equation for the 
density matrix after the Born approximation. In the in¬ 
teraction picture this has the general form: 

kl '' 


where Ok (t) is an operator in the interaction picture and 
VkiiT) accounts for the system-bath coupling, and the 
integral over the bath density of states. From this one 
may derive the population equation 


kl ' 


dt'r]ki{t-t') ( -0|(t)^ (f),6fe(t) 


where (...)/ = Tr[... (<')]. The BR population equa¬ 
tion then follows by assuming has a slow time de¬ 

pendence, and performing the integral over dt' account¬ 
ing only for the time dependence of the interaction pic¬ 
ture operator 

If we focus on late times this procedure is somewhat 
strange, as it is clear that for a problem which has a 
time-independent Hamiltonian in the Schrodinger pic¬ 
ture it is the density matrix in the Schrodinger pic¬ 
ture which will be time independent. As such, an al¬ 
ternate procedure suggests itself: to consider p^F(t') = 




The explicit time dependence of this 
density matrix can be eliminated using Tr[Op*^'^)(f')] = 


Tr 


p(S) 


SO that we have: 


dtF,, = 


kl ' 


drpki (r) ( 'ipl (r)V- ■ (r), Ofc (r) 


.Oi 


where ■ ■)s = Tr[... and we have written t = t — 
t'. Following this prescription, one can again find an 
equation for the vector of real quantities f in the form 
dti — —M'^f + fo, but the matrix M‘® has a different 
form. The matrix is now given by: 


M® = 


f 2 <pIK 

0 

2(fa‘PbK 

\-2ipa<PbK: 


0 

MK'b 

2ipb(PaKb 


^a^bK ^aVbK\ 

‘Pb’PaK'b 


rf 

Ei 


-E^ 

rf 


/ 


( 21 ) 


where now = (wf, - (fllK”) - (uja - and F^ = 

‘PaE'b + ^lE'a- For want of a better name, we refer to 
this as the Schrodinger picture Bloch-Redfield (SpBR) 
equation. The constant vector fo is unchanged. 

The difference between the BR and SpBR equations 
has a simple structure: it corresponds to swapping which 
frequency the bath is to be sampled at in the third and 
fourth column. The origin of this change is the uni¬ 
tary transformation between the interaction and 

Schrodinger pictures, which has the effect of swapping 
time dependence of some “off-diagonal” terms. These 
small changes to the matrix M have several remarkable 
consequences. Firstly we may immediately check that 
the sum rule as written in Eq. (20) is now exactly sat¬ 


isfied. Secondly, one may also consider the behavior of 
the eigenvalues of Eq. (21). Unlike Eq. (13), there is no 


simple closed-form expression for the eigenvalues in the 
general case — Eq. (13) was special in having a structure 


that the secular equation could be written as a quadratic 
in — "^[Ka -I- Kb], but this does not hold for Eq. (21). 
However, one can perform perturbation theory around 
the point A = 0. Clearly the eigenvalues of M and 
M'^ match at this point, as the only distinctions occur if 
Ka 7 ^ Kb- Thus, using standard (non-self-adjoint) per¬ 
turbation theory in terms of the small parameters A and 
Ka — Kb one hnds the lowest SpBR eigenvalue takes the 
form: 


Mo = 




0{A^ 


( 22 ) 


Remarkably, this is identical to the exact solution, fur¬ 
ther confirming the idea that this SpBR equation is an 
improvement over the BR population equations discussed 
previously. 

As we have already seen above, the BR master equa¬ 
tion matches the exact solution well as long as damping is 
weak enough and the Markov approximation is well justi¬ 
fied. The decay rates near resonance have further shown 
that while the BR master equation is correct to leading 
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order in the damping rate, the SpBR equation is cor¬ 
rect to higher order. This suggests that as the damping 
rate becomes larger, the SpBR may give a better numeri¬ 
cal agreement with the exact solution. This is indeed the 
case, and is shown in Fig.j^where we compare the steady 
state values of Fij. Comparing the coherence Fab{t), it 
is clear the SpBR matches the exact solution better than 
the BR approach. The lower panel shows that the two 
theories give very similar results for the populations. For 
the parameters corresponding to Fig.[^the BR and SpBR 
lines would be indistinguishable. 

Note that at A = 0, the SpBR and BR formalisms are 
identical, and so one might expect the results to match 
at this point. However, the matrices M are singular at 
A = 0 (as seen earlier from their eigenvalues). As such, 
the finite population and coherence at A —> 0 correspond 
to a singular limit. 



Q DD _I_ DD I _£__ I DD _ 

-0.1 -0.05 0 0.05 0.1 

-A=(c0b-C0a)/2 


FIG. 3. (Color online) Comparison of steady state values of 
Fij between the exact (solid), BR master equation (dashed) 
and SpBR master equation (dotted), plotted for a larger bath 
density of states Jo = 0.02 and all other parameters as for 
Fig-i 


VI. EXTENSION TO MULTIPLE BATHS 


Extending either the exact solution or the BR mas¬ 
ter equation to multiple baths is simple. For the BR 
master equation, one just finds a separate set of Lamb- 
shift terms h jj a nd dissipator terms LF for each bath, 
so that Eq. (10) involves a summation over contribu¬ 
tions from the baths. Similarly, the expressions for the 
matrix M follow as before, but now with a sum over 
baths, and even the analytic form of the eige nvalues re¬ 
mains true, with Ki i—in Eq. (14). The ex¬ 
act solution is however more complicated. Equation (§ 
still holds, however there is a sum over baths, and each 
term now acquires a bath label: J{v),nB{v), W^(v,t) !—>■ 

The last of these quantities 


now has a more complicated form 


wt\Ft) = 


dC e-<* 

2tt i' — — iO 


T.e.iK'ivf 


(23) 


where the matrix Q can be defined in terms of its inverse, 

m)-% = - c - * 0 ) + 

In the presence of multiple baths, the singular behavior 
at oja = oJb no longer occurs — one may check this by cal¬ 
culating the zeros of Det[0((C)“^]: one now finds there is 
no longer a zero mode, unless the coefficients happen 
to be parallel for different n. The physical origin of this is 
that with multiple linearly independent baths there is no 
longer a linear combination of fields ipi which decouples 
from the baths, and so all modes are damped. As such, 
the collective dephasing model is never correct for pre¬ 
dicting the steady state coherence. The non-secularized 
BR approach continues to correctly describe the system 
as one varies detuning. 


VII. CONCLUSIONS 

In conclusion we have compared the exact and Bloch- 
Redfield solutions for a system of two bosonic modes cou¬ 
pled to a common bath. The late-time behaviors show 
singular dependence on detuning: exactly on resonance, 
significant coherence exists at late times, but for arbi¬ 
trarily small detuning the coherence drops to a smaller 
value which depends on the frequency dependence of the 
density of states. This singular limit appears only at late 
times, corresponding to a slow decay rate for coherence 
that vanishes at the degenerate point. All aspects of this 
behavior are reproduced correctly by a non-secularized 
Bloch-Redfield theory, whereas secularization leads to 
incorrect predictions. The Bloch-Redfield theory does 
not guarantee positivity, nonetheless one can prove that 
the equations describe bounded dynamics of physical ob¬ 
servables, as long as the Markov approximation remains 
valid. A modification to the Bloch-Redfield theory — as¬ 
suming it is the Schrodinger picture density matrix that 
evolves slowly, rather than the interaction picture one 
— leads to an improved time-local theory which satisfies 
required sum rules and exactly matches damping rates 
near resonance. 
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